#include<iostream>
#include<fstream>
extern "C"{
#include<math.h>
#include<time.h>
#include<unistd.h>
#include<sys/time.h>
}
#include "IMParton.h"
using namespace std;

int main()
{
    struct timeval tstart,tend;
    gettimeofday(&tstart,NULL);       //get the start time

    IMParton proton(1,1);             //  Z=1 and A=1, the nuclei is proton
    double lnxstep=(log(1)-log(0.1))/100;       //step length of log(x) for output data
    ofstream outdata1("Q2_0.5_xgluon.dat");    //output file for xgluon at Q^2 = 0.5 GeV^2
    ofstream outdata2("Q2_100_xgluon.dat");    //output file for xgluon at Q^2 = 100 GeV^2
    ofstream outdata3("Q2_0.5_xuv.dat");       //output file for xuv at Q^2 = 0.5 GeV^2
    ofstream outdata4("Q2_100_xuv.dat");       //output file for xuv at Q^2 = 100 GeV^2
    //get the parton distribution functions
    for(int i=0;i<901;i++)outdata1<<exp(log(1e-9)+i*lnxstep)<<" "<<proton.getPDF(0,exp(log(1e-9)+i*lnxstep), 0.5)<<endl;
    for(int i=0;i<901;i++)outdata2<<exp(log(1e-9)+i*lnxstep)<<" "<<proton.getPDF(0,exp(log(1e-9)+i*lnxstep), 100)<<endl;
    for(int i=0;i<901;i++)outdata3<<exp(log(1e-9)+i*lnxstep)<<" "<<proton.getPDF(1,exp(log(1e-9)+i*lnxstep), 0.5)<<endl;
    for(int i=0;i<901;i++)outdata4<<exp(log(1e-9)+i*lnxstep)<<" "<<proton.getPDF(1,exp(log(1e-9)+i*lnxstep), 100)<<endl;

    gettimeofday(&tend,NULL);         //get the end time
    double timeuse=1000000*(tend.tv_sec-tstart.tv_sec)+(tend.tv_usec-tstart.tv_usec);
    //display the total program running time
    cout<<"    Runtime : "<<(timeuse/1000.0)<<" ms."<<endl;
}
